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ABSTRACT 


The  general  purpose  of  the  work  reported  here  is  to  obtain  the  best 
possible  feasible  signal  processing  algorithms  for  estimating  from  radar 
data  characteristics  and  trajectory  parameters  of  bodies  moving  in  the 
air. 


By  analyzing  the  characteristics  of  a  family  of  projectiles,  a  set 
of  approximate  simplified  dynamics  equations  is  obtained  (Section  2) 
that  can  be  used  for  extrapolation  and  backtracking  with  radar  determin - 
able  parameters.  The  formulation  includes  terms  which  account  for  both 
drag  and  drift.  Tentative  numerical  results  indicate  small  resultant 
backtracking  errors  due  to  drift  but  somewhat  larger  errors  due  to  drag, 
especially  when  the  projectile's  velocity  passes  thru  mach  1  during  the 
observation  period. 

In  Section  3>  formulation  is  developed  which  can  be  programmed  to 
determine  the  effects  of  random  radar  errors  on  trajectory  estimation 
and  backtracking  accuracy. 

Section  4  presents  the  formulation  necessary  for  constructing  an 
extended  Kalman  filter  and  smoother  algorithm  for  extracting  pertinent 
state  vectors  and  trajectory  parameters  from  radar  data.  Specific  in¬ 
puts  required  and  relationships  are  given  in  form  suitable  for  programm¬ 
ing. 

Six  appendices  are  incD.uded,  giving  some  of  the  mathematical  de¬ 
tails,  describing  modifications  and  other  experience  with  the  modified 
point-mass  dynamics  program,  and  a  glossary  of  the  main  symbols  used. 
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FOREWORD 


This  report  describes  work  done  from  1  August  1971  through  31  October  1971 
at  the  Moore  School  of  Electrical  Engineering,  University  of  Pennsylvania,  under 
contract  number  DAAB07-71-C-0212  with  U.  S.  Army  Electronics  Command  for  re¬ 
search  entitled  "The  Estimation  of  Trajectories  from  Radar  Data".  The  cognizant 
technical  personnel  at  USA  ECOM  are  Dr.  Leonard  Hatkin,  head  of  the  Radar  Techni¬ 
cal  Area  and  Mr.  Duane  Sheppard,  CSTA  and  31  Laboratory,  Evans  Area  Fort  Monmouth 
N.J.  07703- 
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1.0.  IRTROBUCIipjr.  ‘  *  :  ' 

As  indicated  in  the  previous  Quarterly  Progress  Report,  this  work  is  con¬ 
cerned  with  the  general  signal  processing  problem  of  estimating  characteristics 
and  trajectory  parameters  of  a  body  moving  in  the  air  from  noisy  radar  observa¬ 
tions.  >  !  * 

.  I 

1  The  main  thrust  during  this  last  quarter  has  been  the  adaptation  of  a  BRL 
program  supplied  to'  us  by  USA2C0M:  to  provide  trajectory  data  for  a  suitable 
family  of  ballistic  projectiles,  to  calculate  simulated  radar  coordinate  data 
that  can  result  from  different  geometrical  configurations,  and  to  form  the 
basis  of  "simplified  dynamics "  equations  to  be‘  used  in  filter-smoothing  algorithm 
design  and  for  backtracking  simulation  and  evaluation.  i 

i  ■  > 

Additional  results  have  been  obtained  Concerned  with  evaluation  of  irreduc¬ 
ible  estimation  errors  due  to  radar  measurement  errors  and  the  development  of 
the  necessary  formulation  for  the  optimal  smoothing,  filtering,  and  backtracking 
algorithms  for  the  estimation  of  necessary  trajectory  parameters. 


I 


1. 1  Summary,  main  results. 


Section  2. 0  to  follow  ’describes  the  development  of  simplified-dynamics  equa¬ 
tions,  that  are  suitable  for  approximate  trajectory  calculations  in  terras  of 
.  measurable  parameters.  Starting  with  the  BRL  modified  point-mass  equations  and  a 

1  representative  set  of  projectile  characteristics,  simplified  expressions  for  drag 
•  force  and  for  drift  for’ce  were  evolved.  For  drag  a  universal; coefficient  curve 
was  calculated  that  together  with  a  projectile  variable  scale  parameter  that  must 
i  -  ,  be  estimated  enabled  the  approximate  drag-  force ( to  be  calculated.  For  drift,  a 

normalized  parameter  differential  equation  with  a  universal  curve  coefficient 
was  postulated  requiring  an  initial  ‘value  estimation  from  the  data.  Magnus  force 
}  was  ^assumed  negligible.  In  Section  2,  the  necessary  formulation  is  developed  and 

a  i^et  of  tests  are  described  which  indicate  the  effectiveness  of  the  simplified- 
dynamics  ’equations  in  backtracking  along  the  trajectory  to  the  launch  point.  It 
is  shown  that  drift  force  errors  are  effectively  accounted  for  but  drag  forces 
stiU.  lead  to  errors  that  nre:  perhaps  larger  than  desirable.  Associated  with 
,  this  section  are  Appendix  A,  giving  the  BRL  Modified  point  mass  equations.  Appen¬ 
dix  B  showing  the  Universal  Aerodynamic  Functions  evolved  for  drag  and  drift 
;  force  effects  and  Appendix  C  commenting  on  the  computational  aspects  of  the  work, 
i 


Section  3-0  presents  the  formulation  inecessary  to  determine  theoretical 
limitations ■ on  accuracy  of  trajectory  parameter  estimates  (launch-point  state 
and  aerodynamic  parameters)  due  to  random  rad.x-  noise  errors.  Associated  with 
this  section  is  Appendix  D  which  gives  the  formulation  required  for  transforma¬ 
tions  between  radar -centered  arid  gun-centered  coordinate  systems.  Formulation 
to  account  for  radar  bias  errors  will  be  presented  in  a  subsequent  report. 


■  Section  4.0  presents  the  general  filter- smoother  algorithm  developed  to 
estimate  trajectory  parameters  from  the  noisy  radar  data.  Following  the  ex¬ 
tended  Kalman  filter  technique  an  augmented  time-varying  state-vector  is 


i-l- 


J 


I 


I 


? 


defined  having  eight  components  consisting  of  target  position  ana  velocity  vector 
components  as  well  as  drag  and  drift  parameters.  The  formulation  is  set  up  for 
calculating  optimal  estiroa-o.es  of  present  value  (filter)  and  initial  value 
(smoother)  of  the  augmented  state  vector  updated  on  a  quasi-real  time  basis. 
Initial  work  assumes  availability  of  radar  data  for  off-line  processing  but  once 
algorithms  are  efficiently  programmed,  the  feasibility  of  true  real-time  pro¬ 
cessing  will  be  assessed.  Provision  is  made  for  use  of  simulation  as  well  as 
analytical  techniques  for  performance  and  feasibility  evaluation  of  the  algorithms 
developed.  Results  are  presented  in  a  form  convenient  for  programming.  Specific 
definition  of  matrix  elements  required  is  given  in  Appendix  E.  A  detailed  theore¬ 
tical  development  will  be  presented  in  the  next  quarterly  report. 

Conclusions  and  plans  are  described  in  Section  5. 


2.0  SIMPLIFIED-DYNAMICS  EQUATIONS. 

BRL  (B»Uistics  Research  Laboratory)  Report  No.  1314. ,  Reference  [l],  pre¬ 
sents  the  modified  point-mass  equations  used  in  computing  firing  tables.  These 
equations  are  described  briefly  in  Appendix  A. 

In  order  to  compute  backward  in  time  from  the  observed  projectile  positions 
to  the  launch  point,  some  equations  governing  the  projectile  motion  must  be  used. 
If  the  projectile  type  and  its  aerodynamic  characteristics  were  known,  the  modi¬ 
fied  point-mass  equations  could  be  used  for  computing  backward  along  the  trajec¬ 
tory.  Under  the  conditions  of  interest  in  the  present  study,  nothing  is  known 
about  the  projectile  except  what  is  deduced  from  its  observed  behavior.  It  is 
assumed  in  the  present  study  that  the  accuracy  and  reproducibility  requirements 
on  artillery  fire  must  result  in  generally  similar  aerodynamic  characteristics 
for  all  artillery  shells,  so  that  the  105mm,  155mm,  175mm,  and  6  inch  projec¬ 
tiles  —  called  for  brevity  the  standard  projectiles  —  have  characteristics 
generally  similar  to  any  effective  artillery  shell.  In  particular  it  is  assumed 
that  the  standard  projectiles  can  be  imbedded,  in  a  sense  made  specific  below  in 
connection  with  equations  (5)  and  (12),  in  the  family  of  all  projectiles,  the 
different  members  of  the  family  differing  from  one  another  only  in  numerical 
values  of  certain  parameters. 

Casual  examination  of  the  point-mass  equations  in  Appendix  A  shows  that  the 
trajectory  is  affected  by  three  aerodynamic  forces,  causing  drag,  dr  n't,  and  the 
Magnus  effect,  in  addition  to  gravity  and  the  Coriolis  ''force".  Among  them  they 
involve  seven  functions  of  Mach  number  and  at  least  five  constants,  that  can  vary 
from  projectile  to  projectile.  To  determine  all  these  functions  and  constants 
from  observation  of  the  behavior  of  an  unknown  projectile  would  require  accura¬ 
cies  in  and  durations  of  radar  observations  far  in  excess  of  what  is  achievable. 
The  present  study  assumes  that  it  will  be  enough  to  estimate,  from  the  observed 
projectile  behavior,  the  principal  (zero-yaw)  component  of  the  drag  acceleration 
and  the  drift  acceleration.  The  conditions  under  which  this  assumption  is  valid 
are  discussed  in  2.4  below.  In  all  previous  work  of  which  the  writers  are  aware, 
only  the  zero-yaw  drag  acceleration  has  been  estimated. 


2.1  Drag  acceleration.  The  direction  of  the  drag  force  is  in  the  direction 
opposite  to  that  of  the  instantaneous  velocity.  Its  magnitude  is,  for  the  four 
standard  projectiles, 

(p/ci)KDi(M)v2  i  =  1,2, 3*4  .  (1) 

In  equation  (l)  K-^M)  is  a  function  of  Mach  number  characteristic  of  the  pro¬ 
jectile,  p  is  the  density  of  the  air  mass,  and  v  is  the  projectile  air  speed. 
The  ballistic  coefficient  C,  is  an  empirical  constant  approximately  equal  to 

ratio  m/d  but  generally  also  a  function  of  firing  charge.  Examination  of  the 
four  Kp^(M)  functions,  which  are  illustrated  in  Appendix  B  (actually  the  ratios 

K^/c^  are  plotted  for  nominal  values  of  (I ),  shows  that  the  functions  are 

BRL  data  for  the  standard  projectiles  indicate  that  for  the  155mm,  175cm,  and 
3  inch  shells  is  a  function  of  charge)  is  also  a  function  of  quadrant 

elevation  for  the  155mm  projectile. 
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approximately  proportional,  i.e.  t'he.  ratio  of  any  two  functions,  or  drag  curves, 
is  a  constant  nearly  independent  of  Mach  number.  It  is  this  characteristic  of 
the  functions  that  has  simplified  the  problem  of  estimating  the  drag  effects 

for  an  unknown  projectile.  It  permits  the  selection  of  a  function  Kd(m),  -which 

will  be  devoted  as  the  universal  drag  function >  that  is  independent  of  projec¬ 
tile.  The  universal  curve  has  the  property  that  the  ratio 


(K^CmVc.V^M)  i  =  1,2,3,4 


is  nearly  constant  with  Mach  number,  and  can  therefore  be  approximated  by  a 
constant 


ei =  (KjjiM/c.yyM) .  o) 

Substituting  equation  3  into  equation  (1),  we  have  that  the  drag  acceleration 
for  the  standard  projectiles  is  approximately 

p  c^M)/  (4) 

The  assumption  mentioned  above  that  the  four  standard  missiles  can  be  imbedded 
in  the  family  of  all  projectiles  means,  with  respect  to  zero-yaw  drag  accelera¬ 
tion,  that  for  any  projectile  the  drag  acceleration  will  be  representable  by 

cK^Mjv2  (5) 

where  the  value  of  c  is  to  be  estimated  from  the  observed  projectile  motion. 

The  function  Kjj(M)  is  selected  as  described  in  2.3  below. 

2.2  Drift  acceleration.  It  can  be  seen  from  Appendix  A  that  the  vector  accel¬ 
eration  due  to  drift  is 

-(a^CMjN/v2)^  x  £)  (6) 

for  the  i  standard  projectile,  with  (M)  a  function  of  Mach  number,  N  the 

spin  angular  velocity,  %  the  vector  of  projectile  velocity  relative  to  the 
ground.  As  indicated  in  the  appendix,  the  coefficient  a  is  a  function  of 
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several  physical  parameters  characteristic  of  the  projectile  and  the  so  called 
"lift  factor"  &}  a  scale  factor  function  of  the  firing  charge. 

The  spin  angular  velocity  is  governed  by  the  expression 


N  =  -bipKAi(M)v  H 


(7) 


where  b^  is  again  a  projectile  characteristic  constant  and  K^(M)  another  func¬ 
tion  of  Mach  number.  Evaluation  of  the  constants  K  in  (7)>  using  data  not 

shewn  in  Reference  [l]  but  given  in  ERL  data  associated  with  the  computer  pro¬ 
grams  for  the  modified  point-mass  computations,  shows  that  the  values  for  the 
four  standard  projectiles  are  0.0255#  0.0193#  0.0166,  and  0.0157  ffc2/pound. 
Similar  evaluation  of  the  functions  X^CM)  shows  that  the  same  function  is  given 

by  BRL  for  three  of  the  four  projectiles*  its  value  is  0.007  (dimensionless)  at 
M  =  0  and  it  falls  smoothly  to  about  half  the  value  at  M  =  2.5*  For  the  8-inch 
projectile,  the  value  is  given  as  0.005,  constant  with  Mach  number. 

Clearly  the  situation  is  less  favorable  for  (7)  than  for  the  drag  accelera¬ 
tion,  with  regard  to  embedding  the  four  standard  projectiles  in  some  general 
functional  form  with  small  relative  errors  of  approximation.  Fortunately  the 
drift  effects  are  smaller  than  the  drag  effects  to  start  with,  and  larger  rela¬ 
tive  errors  are  tolerable.  In  the  simplified-dynamics  equations  it  is  assumed 
that 


N  =  -0.020pKA(M)vN  (8) 

is  a  tolerable  universal  approximation  to  (7)#  with  Ka(M)  taken  as  the  function 
K^M)  for  the  105mm,  155mm,  and  YlJm  projectiles. 

The  expression  (6)  may  be  written 

-  ( M )  N/v2 )  (x  *  S) 

where,  as  explained  in  2.3  below,  a  function  K(M)  and  four  constants  are  used 
(these  c^’s  are  unrelated  to  the  c^s  of  (3)  and  (4)*  they  merely  serve  an  ana- 


i  =  1,2, 3, 4 


(9) 


logo  us  purpose)  such  that  for  eachT”i  -  1,2, 3, 4,  K(24)  approximates  c .  K.  (M)  as 

closely  as  possible.  Again,  the  relative  accuracy  will  be  poor  by  comparison 
with  the  relative  accuracies  achieved  for  the  drag  approximations,  but  as  will 
be  seen  in  2.4  below  the  use  of  the  drift  terms  in  the  simplified-dynamics  equa¬ 
tions  gives  results  considerably  better  than  does  ignoring  the  drift. 

The  drift  acceleration  can  now  be  taken  as 

..(ai/ci)K(M)(N/v2)(2  x  £)  i  =  1,2. 3,4  (10) 


where  a^/c.  is  a  single  parameter  and  ic(M)  a  universal  function,  and  where  U  is 

given  by  (§).  As  in  the  passage  from  (4)  to  (5),  we  could  imbed  the  values  of 
the  parameter 


r 


i 


(11) 


in  a  family  of  parameter  values,  and  for  the  case  of  an  unknown  projectile  esti¬ 
mate  the  value  of  r  from  observed  drift  acceleration  effects  assumed  to  be 
governed  by  an  expression 


-rK(M)(K/v2)(y  X  i) 


(12) 


together  with  (8),  in  which  only  r  needs  to  be  estimated.  However,  although 
(8)  through  (12)  are  useful  for  exhibiting  the  underlying  reasoning,  it  is 
more  convenient  in  the  computations  to  define  a  new  variable  s(t)  to  be  estimat¬ 
ed,  as  part  of  the  state  vector,  than  to  estimate  the  parameter  r. 

Let 


s(t)  =  rN(t)  =  6N(t)/c 


(13) 


by  definition.  Then  (8)  becomes 


s  =  -0. 020pK. (M)vs 


(14 


and  the  drift  acceleration  (12)  becomes 


-(s/v2)K(M)fe  X  i) 


(15) 


Equation  (14)  and  expression  (15)  in  the  simplified-dynamics  equations  replace 
(7)  and  (6)  of  the  modified  point-mass  equations. 


^•3  Choice  of  universal  functions .  In  passing  from. (9)  to  (lo),  it  is  necessary 
to  find  a  function  K(M)  and  four  constants  c±,  i  =  1,2,3,4,  so  that  each  of  four 

given  functions  K±(M)  is  approximated  by  c±K(m).  Below  is  described  the  basis 
for  choosing  the  universal  function  ic(M)  and  the  four  constants  c  .  The  same 

problem  occurs  in  passing  from  (3)  to  (4),  with  C.K_,  (m)  given  instead  of  K4  (M). 
i  =  1, 2,3,4.  1  m  i 

One  way  to  achieve  the  desired  result  is  to  minimize 

“b 

s  =  ?  jM  wi(M){ciK(M)  -  K^Mj^dM  (16) 

Or 


where  (M&,Mb)  is  the  range  of  interest  of  M  and  w^M)  is  any  set  of  weighting 

functions  that  permit  assigning  importance  to  certain  i,  to  certain  subranges 
of  M,  or  both.  In  what  follows,  all  integrals  are  over  (M  ,M  )  and  all  sums 
are  over  i.  a  “ 

For  a  given  set  (c1,...,c^),  the  K(M)  that  minimizes  (16)  must  satisfy 


S(K(M)  +  6k(M))=  S(K(M)) 


(17) 


for  an  arbitrary  small  variation  6k(m),  by  the  reasoning  basic  to  the  calculus 
of  variations.  To  first-order  effects,  (17)  implies 


j2{wi(M)(ci6K(K))(ciic(M)  -  K.(M))  dM  =  0  (lB) 


which  can  hold  for  arbitrary  small  6k(m)  only  if 
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(19) 


2ciwi(M)(c.K(M)  -  K±(M))  =  0 


or 


K(m)  =  2  (M)^  (M)  (M) 


(20) 


Or  if  X(M)  is  fixed,  the  c’s  that  minimize  must  satisfy  Bs/oc.  =  0, 
i  =  1,2, 3,4,  or  x 

2j‘wi(M)K(M)(c1K(M)  -  K^MjJdM  =  0  i  =  1,2, 3,4  (21) 

from  which 

C±  «  Jwi(M)K(M)Ki(M)dM/jw1(M)(K(M))2dM  i  =  1,2, 3,4  (22) 


By  applying  (20)  and  (22)  alternately,  since  each  application  reduces  S, 
one  can  find  a  succession  of  c^’s  and  functions  K(M)  until  successive  changes 

in  S  are  as  small  as  desired.  However  there  is  no  unique  solution  to  theT  problem 
as  stated  above,  because  for  any  real  k^O,  multiplication  of  each  c  by  k  and 
division  of  k(M)  by  k  leaves  S  unchanged.  It  is  convenient  to  make  the  solution 
unique.  This  is  done  by  using  the  following  computation  procedure: 

(a)  Take  an  arbitrary  f(M)  to  begin  with  —  for  example  K^M). 

(b)  Use  (22),  with  f(M)  in  place  of  K(M),  to  compute  a  set  of  c’s. 

(c)  Normalize  the  c’s  by  dividing  each  by  their  mean,  so  that  the  normalized  c’s 
will  satisfy  2c^  =»  4. 

(d)  Using  the  normalized  c’s,  vise  (20)  to  compute  ic(M). 

(e)  Use  (22),  and  normalize  as  in  (c),  to  compute  a  new  set  of  c’s. 

(f)  Repeat  (d)  and  (e)  until  the  changes  in  an  iteration  are  satisfactorily  small. 
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She  above  procedure  was  used  with  the  w^(M)  =  1  for  all  M  and  i  in  order  to 
determine  the  universal  drag  function  K^(M)  and  the  corresponding  four  constants 
ci*  The  nominal  ballistic  coefficients  for  the  155mm,  175mm,  and  8  inch  projec¬ 
tiles  and  the  correct  single  value  for  the  1052m  projectile  given  in  Appendix  B 
were  assumed  for  the  computation.  It  should  be  noted  that  the  above  choice  of 
weighting  functions  although  an  obvious  one,  is  quite  arbitrary.  Other  choices 
could  have  been  made,  possibly  in  order  to  shape  the  universal  curve  to  minimize 
drag  computation  errors  in  the  transonic  region  for  specific  expected  target  pro¬ 
jectiles. 

Only_one  cycle  of  the  procedure  was  used  for  determining  the  universal  drift 
function  ic(M)  and  the  associated  constants  r^.  Justified  by  our  assumption  that 

large  relative  inaccuracies  are  acceptable  in  drift,  w. (M)  was  taken  to  be  the 

i 

Dirac  delta  function  6(m-Mq)  with  MQ  =  0.8  for  all  i.  As  was  previously  indicated, 
the  computations  were  scaled  by  the  nominal  lift  coefficients  of  Appendix  B. 


2.4  Test  of  the  simplified-dynamics  equations.  One  limiting  factor  in  the  per¬ 
formance  of  the  system  studied  by  Project  RAIRAM  is  the  amoung  of  error  incurred 
in  using  the  simplified-dynamics  equations  for  the  backward  integration  along 
the  trajectory.  To  measure  the  amount  of  this  error,  it  is  intended  to  use  the 
modified  point-mass  equations  to  compute  a  trajectory  up  to  some  point  (x,y,z,t), 
and  then  to  use  the  simplified-dynamics  equations  for  computing  backward,  as 
these  equations  will  be  used  in  practice  for  computing  the  estimated  launch 
point.  The  difference  in  x  and  z  coordinates  between  the  point  on  the  backward 
trajectory  at  which  y  =  0  and  the  starting  point  for  the  forward  trajectory  com¬ 
putation  is  the  error  due  to  the  use  of  the  simplified-dynamics  equations,  plus 
numerical  errors  in  computation;  the  effects  of  the  latter  are  estimated  by  an 
auxiliary  computation  in  which  the  modified  point-mass  equations  are  used  in  the 
backward  computation  as  well  as  in  the  forward  computation. 


In  the  backward  computation,  the  correct  parameter  value  is  used  in  the 

drag  computations  and  the  correct  starting  value  is  used  for  s(t),  the  drift 
variable.  In  both  cases  the  constants  were  properly  scaled  to  account  for  the 
differences  between  the  nominal  drag  and  lift  coefficients  of  Appendix  B  and 
true  values  for  the  particular  projectile  cases  considered.  There  will  be  errors, 
in  practice,  in  the  estimates  of  c.  and  s(t),  but  these  errors  are  not  chargeable 
to  use  of  the  simplified  dynamics  equations. 

The  forward-backward  computation  test  on  the  simplified-dynamics  equations 
has  been  performed  for  several  trajectories;  results  are  shown  in  Table  1.  Char¬ 
acteristics  of  the  test  trajectories  used  are  shown  in  Figure  1. 
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TABLE  1  BACKTRACKING  ERRORS  FOR  A  SELECTION  OF  TRAJECTORIES 
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ERRORS  IN  BACKTRACKING  TO  LAUNCH  POINT 
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FIGURE  1  Characteristics  of  Test  Trajectories 


The  table  summarizes  the  results  of  three  kinds  of  computations.  In  all 
three,  the  BRL  modified  point-mass  program  was  used  to  compute  trajectories  to 
some  preselected  time,  the  ’time  at  reversal”.  From  the  state  parameters  at 
the  time  of  reversal,  three  backtracking  computations  were  made,  proceeding 
backward  in  time  to  zero  elevation.  The  three  differ  as  follows.  In  the  case 
of  most  interest,  the  backward  computation  used  the  simplified  dynamics  equa¬ 
tions;  the  results  are  given  in  the  columns  headed  "simplified  dynamics  includ¬ 
ing  drag  and  drift".  Magnus  effects  were  not  included  in  the  backward  computa¬ 
tions.  In  order  to  show  what  would  happen  if  no  drift  effects  are  included, 
drift  force  was  set  equal  to  zero  in  another  set  of  backward  computations, 
reported  under  "simplified  dynamics  including  drag  only".  Finally,  to  estimate 
errors  due  to  numerical  round-off,  a  third  set  was  made  using  the  complete  BRL 
equations;  the  results  are  shown  under  "BRL  equations". 

The  tables  show  that,  for  the  trajectories  tried,  (a)  the  effects  of  com¬ 
putational  errors  are  negligible,  (b)  drift  is  the  principal  source  of  error  if 
only  drag  effects  are  included  in  the  backward  computations,  and  (c)  drift  error 
effects  can  be  reduced  to  amounts  small  by  comparison  with  drag  error  effects, 
insofar  as  errors  in  the  simplified-dynamics  equations  are  concerned.  (Whether 
the  filter  can  estimate  the  drift-variable  well  enough  to  take  advantage  of  the 
situation  remains  to  be  seen). 

For  the  trajectories  so  far  tested,  the  errors  due  to  the  use  of  the  simpli¬ 
fied  dynamics  equations  are  below  25  meters  in  estimated  launch  point  for  all 
trajectories  in  which  the  radar  observations  can  be  made  within  20  second  of 
firing. 


3-  0  IRREDUCIBLE  LAUNCH-POINT  ESTIMATION  ERRORS  DUE  TO  RANDOM  RADAR  ERRORS. 


3.1  Purpose.  On  the  assumption  that  the  simplified-dynamics  equations  represent 
the  true  behavior  of  the  projectile,  this  section  describes  how  to  compute  the  co- 
variance  matrix  of  the  error  in  launch-point  estimation  for  least  squares  estima¬ 
tors.  The  simplified-dynamics  equations  are  defined,  ana  the  errors  incurred  in 
their  use  discussed,  in  section  2.  For  the  purposes  of  this  section,  it  is  enough 
to  say  that  they  involve  a  drag  variable  and  a  drift  parameter  to  be  estimated. 

It  is  assumed  here  that  the  only  source  of  error  is  the  random  radar  errors. 
Bias  radar  errors  will  be  investigated  separately.  The  errors  in  launch-point 
estimation  here  studied  will  depend  only  on  the  covariance  matrix  of  the  radar 
errors,  the  trajectory-radar  geometry,  and  the  schedule  of  radar  observations. 

They  are  sampling  errors,  in  the  sense  that  if  enough  independent  radar  observa¬ 
tions  could  be  made,  the  errors  in  launch-point  estimation  due  to  the  radar 
errors  could  he  made  as  small  as  desired. 

This  section,  then,  tells  how  to  compute  the  performance  of  the  best  possible 
(least-squares)  signal-processing  algorithms  in  coping  with  random  radar  errors, 
in  deducing  what  needs  to  be  known  about  the  projectile  and  its  trajectory.  Later, 
when  the  performance  of  individual  algorithms  (e.g. ,  Kalman  filters)  is  investi¬ 
gated,  the  results  of  the  computations  described  below  will  provide  limits  on  the 
performance  of  the  algorithms.  The  limits  could  be  reached  if  the  covariance 
matrix  of  the  radar  errors  is  known  and  if  enough  radar  observations  could  be  made 
and  if  computation  time  and  computer  memory  were  not  limited,  comparison  between 
the  limiting  performance  here  considered  and  actual  algorithm  performance  will  be 
used  in  guiding  decisions  on  whether  to  accept  performance  of  a  particular  algorithm 
or  to  seek  improvement  at  the  cost  of  computer  capacity. 


3.2  Outline  of  method.  The  following  is  well-known:  If  there  are  N  observed 
column- vectoru  (radar' observations )  l 


Ik  =  2*  +  k  =  1,...,N  (23) 


where  is  a  true  value  and  each  6r  is  a  zero-mean  normally- distributed  error 

vector  independent  of  the  others,  and  if  there  is  to  be  estimated  a  vector  v 
(augmented  initial  state  vector)  that  satisfies 


2k  "V 


kvi,...,N  (24) 


for  given  matrices  A,  ,  then  the  maximum- likelihood  estimator  for  v  is 

'  5  — 


*  „I  ,T  „I  ~ 

s  =  s 


-k  :  i 


where  the  superscripts  I  and  T  denote  respectively  the  inverse  and  -transpose 
operations, 


n  m  T 

s  v 


,(?6) 


*k  * Ef6^  62* 1 


is  the  covariance  matrix  for  the  errors  in  r,  . 

-k 

Under  the  assumption  that  the  errors  6r.  ,  and  the  errors  6v  that  they  cause 

in  v,  are  small,  the  RATRAN  launch  point  estimation  problem  can  be  identified 
with  the  linear  problem  just  described  as  follows.  The  vector  v  has  eight  com¬ 
ponents:  the  coordinates  in  gun-axes  xq,  y  ,  and  zq  =  0  and  the  velocity  com¬ 
ponents  xq,  yQ,  and  zq  at  the  time  tQ  of  launch)  a  constant  parameter  c  character¬ 
istic  of  the  projectile)  and  a  drift  parameter  CQ>  the  initial  value  of,  a  variable 

used  in  the  simplified-dynamics  equations  for  drift  computations.  For  notational 
convenience  these  eight  numbers  are  denoted  also  by  v^,...,Vg  respectively.  '  The 

vector  r.  of  radar  observations  at  time  t.  has  either  three  or  four  components, 

" K  4  t  s 

depending  on  whether  doppler  is  not  or  is  used,  the  other  three  being  the  range 
and  two  angles.  (Four  different  radar  types  will  be  considered:  range-a-p  and, 
range-azimuth-elevation,  each  with  and  without  doppler. )  Then  : 


\  -  (*uw) 


where 


a.  <*>.&•  «/*., 

ij  i  '  .  a* 


i 


i 


1 


the.  change  in  the!  ith  radar  coordinate  at  time!  t.  due  to  small  .unit  change  in  the 
jth  component  of  the  initial  augmented  state  vector  v. 

5  !  “  * 

Equation  (25 )‘  implies  that  :  , 


E[6v  6y^}  =  S1 


(30) 


:  l 


is  the  covariance  matrix  for  the  error 


S'  A 

I  =  Z  "  ~ 


(3D 


•  I  A  I  *  * 

in  the  maximum- likelihood  estimator  v.  i  By  the  assumptions  of  normality  and 
linearity,  v  is  al$o  the  least- squares  estimator. 

i  Sihce'the  purpose  of  RATRAN  is  not  to  estimate  v  at  a  particular  value  t  =  0 
of  time,  with  respect  to  the  radar1  observation  timesj  but  rather 'at  a  particular 
i'  value  z  =  0  of  projectile  altitude  —  and  in  general  the1  two  conditions  will  not 
concur  —  v  is_  adjusted  by  . 


or  1 


A  A  A  «  /* 

XT  =  X  -  Z  X  /z 

L  o  ;  o  o'  o 

A  A  A  •  /• 

V.  =  V  -  Z  v  /  z 

o  “o'  o 


,X1 


A 


=  Lv 


(30 


(33) 


with 


r  b 


i 


L  = 


!l  -*  /y  ,  0  0  0,0  0  0 


0"0 

! 


V°  -4A  ’1  0  ,0  0  0  0, 


-15- 


I 


(34) 


i 


Because  the  adjustment  (32),  is  small,  the  true  value  of  L  may  be  used,  and  for 
convenience  will  be.  used,  in  the  computations  here  described,  in  place  of  the 
estimated  value  that  would  have  to  be  used  in  the  actual  radar  system. 

The  covariance  matrix  for  the  error  launch-point  estimate,  using  least- 
squares  estimation,  is  therefore 


=  E(L  v  f  LT)  •■=  LS1  LT 


(35) 


The  computation  of  P  is  the  subject  of  the  remainder  of  this  section. 


3*3  Inputs  to  the  computations.  Each  computation  of  P  is  for  the  following  in¬ 
puts:  A  projectile  trajectory,  requiring  specification  of  the  projectile  and.  its 
initial  velocity  vector)  the  position  of  the  gun  and  its  direction  of  fire  rela¬ 
tive  to  the  radar)  the  choices  with  regard  to  the  radar  —  namely  the  use  of 
doppler,  the  choice  of  angular  coordinates,  the  radar  altitude,  and  (only  if  range- 
a-R  is  used)  the  antenna  array  tilt  angle)  the  covariance  matrices  R.  for  the 

ranuom  radar  errors:  and  the  set  of  observation  times  t  . 

k 

The  projectile  trajectories  will  be  computed  using  the  BRL  modified  point- 
mass  computer  program  as  adapted  for  the  simplified-dynamics  equations.  The  drag 
parameter  v^  =  c  and  the  initial  value  Vg  =  of  the  drift  variable  must  be 

specified.  The  gun  coordinates  z  ,  x  relative  to  the  radar  (y  is  vertical  and  x 

forward,  as  explained  in  Appendix  D),  and  the  azimuth  firing  angle  A2  in  the 
rectilinear  radar-axes  system  must  be  given)  these  numbers  are  specified  as  though 
the  earth  were  flat,  but  the  trajectory  and  the  radar  coordinates  are  computed  for 
a  spherical  earth.  The  covariance  matrices  will  be  specified  either  as  constant 

or  as  a  function  of  radar  coordinates,  never  as  time  functions  as  might  have  been 
inferred  from  the  definition  relative  to  the  observation  times  t^.  The  computer 
program  will  provide  for  the  possibility  that  the  elements  are  to  be  computed 
from  the  radar  coordinates  at  the  time  t^.  The  matrices  R^  will  be  4-by-4,  if 

doppler  is  used)  if  doppler  is  not  used,  the  computations  are  simplified  as  dis¬ 
cussed  in  3-5-3  below.  The  choice  of  observation  times  t^  that  the  program  will 
be  able  to  accomodate  is  discussed  in 


3.4  Values  for  radar  error  covariances.  Inputs  to  this  analysis  program  re¬ 
presenting  the  variances  of  the  random  radar  estimation  errors  can  be  described 
as  follows.  Logical  system  operation  will  consist  of  a  sequence  of  independent 
"looks"  at  each  target,  each  look  consisting  of  a  burst  of  radar  pulses  during 
which  range,  two  angles,  and  range  rate  may  be  estimated.  Aside  from  bias  errors, 
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which  will  be  handled  separately  in  this  study,  the  different  measurement  errors 
can  be  considered  to  be  independent  random  variables  with  zero  mean  values  being 
effectively  independent  also  from  look  to  look.  The  following  formulas  will  be 
applicable  to  this  end  (see  for  example  Skolnik  |_2],  Chapter  10). 


r  2  _  C2  +  Q  2 
R  4vb2(2E/No)  rc 


Range  error  variance 


(36) 


(2E/N0) 


+  o.s 


Angle  error  variance 


(37) 


x2 


4yt2(2E/Nq) 


Range  rate  error  variance  (38) 


In  the  above  expressions,  2E/Nq  is  effectively  signal-to-noise  ratio  deter¬ 
mined  by  transmitted  power,  antenna  gain  and  target  range  and  cress  section;  X  is 
the  wavelength  of  the  transmitted  signal;  C  is  the  velocity  of  light;  y_  is  the 
effective  bandwidth  of  the  transmitted  signal;  \'A  is  the  effective  aperture 
dimension  of  the  radar  antenna;  yt  is  the  effective  duration  of  the  coherent 

™  2  2  2 

signal  used  for  doppler  estimation;  0 _  ,  an  ,  a*  are  the  non-range-dependent 

Rc  °c  Rc 

components  of  measurement  error  variances,  due  to  (unavoidable)  system  imperfec¬ 
tions. 


A  convenient  simplification  is  to  make  a  preliminary  set  of  calculations 
with  a  nominal  system  at  a  nominal  range  R  ,  then  make  use  of  the  fact  that 
2E/Nq  varies  inversely  as  range  to  the  fourth  power  so  that 


2E/No  =  (2E/Nq)o  (R0/R) 


4 


(39) 


Then  equations  (36),  (37)*  &ad  (38)  simplify  as  follows: 
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m 


*  a/(s/R0) 

0 

=  °ea(s/a0>4 

o 

a^(R/Ro)4 


(41) 


2  2  2 

where  aR  ,  cQ  and  a-  are  variances  of  errors  due  to  thermal  noise,  calculated  at 
o  o  Ac 

the  nominal  range  Rq. 

The  above  quantities  are  then  used  as  the  main-diagonal  elements  of  the 
matrix  R^  referred  to  in  section  3*3  above;  all  off-diagonal  elements  being 
zero. 


3-5  Details  of  the  computations.  Let  k  be  arbitrary  but  fixed,  so  that  it  need 
not  be  exhibited  in  the  notation.  Then  (28)  and  (29)  may  he  expressed,  using  nota¬ 
te  on  common  in  control  theory,  by 


T  T 

A  =  dr  /ov.  (4g) 


=  2^  =  ^e  vector  whose  elements  are  the  values  at  time  t^  of  the 

projectile  coordinates  x^,  y^,  z^  and  their  time-derivatives  in  the  gun-axes 
sysls&n* 


2k  "  (xk  yk  zk  \  \  (**) 

Clearly  u  is  determined  by  v,  for  given  t^,  and  r  is  determined  by  u,  and  so 

8r^/ov  =  (8uT/8v)(8rT/du) 
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(45) 


or  A  =  CB  where 


BT  =  duT/dv,  CT  =  dr'Vdu 


(b6) 


in  the  notation  of  (43)* 

The  6-by-8  matrix  B  is  independent  of  the  type  or  location  of  the  radar,  ex¬ 
pressing  only  the  relation  between  the  initial  augmented  state  vector  of  the  pro¬ 
jectile  and  a  later  state  vector.  The  elements  of  B  will  be  computed  and  stored. 
The  computation  of  the  elements  b_  of  B  require  nine  executions  of  the  program 

for  computing  trajectories.  One  of  the  nine  is  called  the  reference  trajectory. 
The  reference  trajectory  begins  with  the  nominal  initial  conditions.  In  each  of 
the  other  eight,  one  of  the  eight  components  v. ,. . . ,Vg  of  the  initial  augmented 

state  vector  is  perturbed 


v.  _  v.  +  6v  j  =  1,...,8  (47) 


For  each  k,  the  values  of 


=  VV 


X  —  1, ... fS 


are  given  by 


=  (l/6^)  (pertui-bed  ^ 


reference  u^) 


m 


with  the  right-hand  side  evaluated  at  time  t^.  Thu  perturbations  6v^  are  input 

constants,  to  be  determined  by  preliminary  trials  not  defined  here  in  detail, 
their  purpose  being  to  determine  values  6v^  for  which  the  differences  between 

perturbed  and  unperturbed  u-values  are  above  the  round-off  error  noise  level, 
but  are  small  enough  for  (48)  tc  approximate  a  partial  derivative.  Forty- 
wight  b’s  will  be  computed  and  stored  for  each  k. 


t 


-19- 


In  addition  to  the  48  matrix  elements  b. there  will  also  be  stored  the  six 
vector  elements  u(t^)  of  (44),  as  computed  fo$  the  reference  trajectory,  for  vise 

in  computing  the  elements 

^  k  -  1,* •  ♦  (49) 

of  Ck>  as  explained  below.  Each  trajectory  thus  requires  storing  54ft  numbers. 

The  matrices  will  not  be  stored  between  runs.  They  will  be  computed  as 
needed,  from  input  constants  and  the  values  of  the  vectors  u(t.  )  =  u^,.  The 

first  step  in  the  computation  of  each  C.  is  the  computation  of  the  projectile 
position  vector 


w,  = 

-k 

yrk 

1 

a  M 

y(\) 

1 

VW 

K/ 

*  n 


(50) 


in  a  right-handed  rectilinear  coordinate  system  with  origin  at  the  radar.  In  (50) , 
the  radar  coordinate  system  has  the  y-axis  along  the  vertical  at  the  radar,  posi¬ 
tive  upward,  and  the  x-axis  along  the  nominal  forward  (zero-azimuth)  radar  direc¬ 
tion.  The  values  of  x(t.  ),  y(t  )  and  z(t  ),  three  of  the  six  elements  of  the 

previously-stored  vector  u^,  are  the  components  of  the  projectile  position  vector 

in  the  right-handed  rectilinear  system  that  the  trajectory  computations  use.  This 
system  has  its  origin  at  the  gun  and  its  y-axis  positive  upward  along  the  local 
vertical,  and  its  xy-plane  contains  the  initial- velocity  vector.  M  is  the  3-by-3 
matrix  of  the  rotation  needed  to  make  the  gun  and  radar  axes  parallel,  and  d  is  the 
position  vector  of  the  gun  in  the  rectilinear  radar  axes  system.  M  and  d  are  com¬ 
puted  in  terms  of  input  constants  as  follows: 


M  =  (m^  )#  =  d.,2, 3 

dT  .  (dx  d2  d3) 

2  2 

s  ■  (2g  +  xg  > 


(51) 

(52) 


(53) 


*See  Appendix  D  for  a  simplified  discussion  of  the  coordinate  conversion  used. 

-20- 


C3  =  -Xg/g  (5*0 

53  =  -Zg/g  (55) 

Cl  =  (S3)  sinA2  +  (C3)  eosA2  (5$) 

SI  =  -(S3)  cosA2+  (C3)  sinA2  (57) 

54  =  g/R,  (58) 

where  R  is  the  radius  of  the  earth. 

c4  =  1  -  (l/2)(s4)2  (59) 

m;L1  =  (C1)(C3)(C4)  -  (SI) (S3)  (60) 

=  -  (C3)(S4)  (61) 

m13  =  -  (Sl)(C3)(C4)  -  (Cl) (S3)  (62) 

=  (C1)(S4)  (ft) 

n>22  =  C4  (6^) 

n^3  =  -  (S1)(S4)  (65) 

m31  =  (C1)(S3)(C4)  +  (S1)(C3)  (66) 

1^2  =  -  (S3)(s4)  (67) 

^  =  -  (Sl)(S3)(c4)  +  (C1)(C3)  (68) 
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(70) 


dg  =  -g(s4)/2  -  hr 

%  =  zg  (71) 

In  (71),  is  the  altitude  of  the  radar  (in  excess  of  the  altitude  of  the  gun). 

The  matrix  M,  with  components  given  by  (6o)  through  (68),  is  also  used  to 
compute 


M 

*rk 

=  M 

y(\) 

\W 

(72) 


The  vectors  and  are  used  to  compute  the  radar  range  at  time  t^ 


r 


k 


(73) 


the  direction  cosines 


the  range  rate 


C5  =  xrk/rk 


06  =  yrk/rk 


C7  =  srk/rk  , 


(C5)i.k  +  (c6)yrk  +  (°7)zrk  > 


and  the  ratios 


(74) 

(75) 

(76) 


(77) 
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CS  =  irk/rk  (T8) 

09  =  yrk/rk  (79) 

CIO  =  zrk/rk  (So) 

By  definition  (see  figs.  1  and  2) 

a  -  arc  gos{(c5)cosy  +  (c6)siny}  (8la) 

where  Y  is  the  antenna  array  tilt  angle. 

p  =  arc  cos(C7)  (82a^ 


The  purpose  of  (50)  through  (82a)  is  the  computation  of  c^^(t^),  i  =  1,...,4- 
j  =  l,...,6j  k  =  1,... ,N,  where  c^  is  the  partial  derivative  of  the  ith  radar  co¬ 
ordinate  with  respect  to  the  jth  component  of  u(t^)  as  defined  by  (44)*  If  the 

first  element  of  r  is  taken  to  be  the  range  r  and  the  fourth  to  be  the  range  rate 
r,  in  both  the  range-a-p  case  and  the  range-azimuth-elevation  case,  then  the  first 
and  fourth  rows  of  C  will  be  the  same  in  the  two  cases.  They  are  given  by 


°ll  '  C5 

(83) 

C32  -  <* 

(84) 

c13  =  07 

(85) 

clk  *  c15  -  C16  “  0 

(86) 

c4l  =  %£  '  °43  ’  0 

(87) 

044  =  08 

(88) 

°45  '  09 

(89) 
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In  the  range -a-p  case. 


Equations  (8la)  and  (82a)'  give  a  and  Performing  the  partial  differentiations 
with  respect  to  the  u^’s  gives 


C21  =  C-((C6)2  +  (c7)2)cosy  +  (C5)(c6)sinY}/DEK0M  (9 

C22  =  UC5)(c6)cosy  -  ((C5)2  +  (C7)2)sinY}/DENOM  (9; 

C23  =  (C7)((C5)cosy  +  (C6)siny)/DEN0M  (9 

C24  =  C25  =  C26  =  0  ^ 

where 

2  2  V2 

DENOM  =  { (xrkeinY  -  y^cosY )  +  zyk  }  (9 

C31  =  (C5)(C7)/RH0  (9 

(9 


C,p  *  (C6)(C7)/RH0 


In  the  range -azimuth-elevation  case. 


2*  -  <rk  “k  s\  V 

(91b) 

where 

* 

tak  =  arc  tan  (trk/xrk) 

(3ib) 

SI^  =  arc  sin  (yrk/rJ 

(32b) 

Performing  the  partial 

differentiations  gives 

//  2  2 
C.,  S'Z./X.  +  Z  . 

2^.  rk'  '  rk  rk 

) 

(92b) 

0 

11 

C\J 

(93b) 

C23  ~  xrk/'xrk  +  2rk  ^ 

(94b) 

C24  =  C25  ”  C2 6  =  0 

(95) 

C31  =  ~Xrk  yrk/rk 

(97b) 

c32  = 

(98b) 

C33  =  "yrk  Zrk/rk^Xrk 

+  2rk2) 

(99b) 

C34  =  C35  =  C36  =  ° 

(100) 

3.5*1  The  observation  times  t^.  As  has  been  said,  there  will  be  first  computed 


and  stored  43  elements  of  and  six  of  u,_  for  time  t 


4 


,t  .  The  matrix 

'  Yt 


n 


T  T 

s  -  kii  K  4  \ 


(26) 


is  approximated,  for  evenly-spaced 

tk  =  \  +  (k-l)f  , 


(102) 
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1 


by 


kJ  Aj/2  +  (lA)  JA  AT  a1  A  dt  +  A  R*  An/2; 


(103) 


where  the  integral  has  been  equated  to  its  trapezoidal-sum  approximation.  Equa¬ 
tion  (103 )  shows  that  the  number  n  of  the  observations  need  net  be  carefully 
selected,  because  for  given  end-points  t^  and  t^, (103)  can  be  used  over  a  range 

of  intermediate  spacirigs  t,  as  long  as  they  are  small  enough  to  warrant  the 
integral-sum  approximation.  Since  there  is  no  reason  to  expect  rapid  changes  in 
the  product  Aa (t)RI(t)A(i)>  it  seems  safe  to  take  the  values  of  t  perhaps  l/2 
second  apart,  for  the  purpose  of  storing  the  elements  of  and  u^,  and  to  use 

the  stored  values  for  computing  the  integral  in  ( 103 )  >  after  which  (103)  may  be 
vised  with  arbitrary  t  to  compute  S  for  other  observation  rates,  by  means  of 


S(t)  =  A?  rJ  A./2  +  AT  R1  A  /2 
'  ill'  n  n  n' 

+  (t/t){s(to)  -  A*  R*  A1/2  -  A*  R l  Aj2]  (104) 


where  t  is  the  spacing  of  interest  and  is  the  spacing  used  in  the  stored  values, 
so  that  S(tq)  is  given  by (26).  Of  course  the  end-points  t^  and  t^  must  be  the  same 
for  the  two  spacings  and  T. 

The  values  of  t^  and  t  will  be  selected  after  some  preliminary  investigation, 

for  each  trajectory,  on  the  earliest  time  the  projectile  will  be  visible  to  the 
radar  and  the  latest  time  for  which  the  accumulated  errors  in  the  simplified- 
dynamics  equations  are  of  tolerable  magnitudes. 

3.5*2  Inversion  of  the  matrix  S.  From  (34)  and  (35)  it  is  clear  that  only  the 
first  three  rows  and  columns  of  S1  are  of  interest.  If  the  8-by-8  matrix  S  is 
partitioned 


k 


S  = 


(105) 
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is  5-by-5,  its  inverse  has  in  its  upper- 


where  S1  is  3-by-3,  Sg  is  3-by-5,  and 
left-hand  corner  the  matrix 


I  T' 
S2  S3  S2, 


)  * 


(106) 


The  computation  of  (lCo)  requires  a  5 -by- 5  inversion  and  a  3-by-3  inversion,  in 
place  of  an  8-by-8  inversion  for  the  entire  matrix  S.  The  use  of  (I06)  is  pre¬ 
ferable,  Then 


replaces  (35 ),  with 


P  =  L  Sx  L 


/I 


-*A 

-*A 


(107) 


doe) 


3.5.3  Effect  of  not  using  doppler.  If  doppler  is  not  used,  the  computations  are 
considerably  simplified.  The  matrices  have  the  following  dimensions: 


Matrix 
S,  S1 
A 

R,  R1 
B 
C 


With  doppler 
rows  columns 

8  •  8 
4  8 

4  4 

6  8 

4  6 


Without  doppler 


rows 

8 

3 

3 

3 

3 


columns 

8 

8 

3 

8 

3 


In  particular,  only  the  first  three  rows  of  B  need  be  computed  and  stored, 
if  it  is  known  that  doppler  will  not  be  used.  It  is  intended  to  use  the  full 
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T 


system,  as  described  above  for  the  case  in  which  doppler  is  used,  at  least  ;for  the 
first  few  trajectories  studied,  to  assess  roughly  the’  value  of  using  doppler.  If 
doppler  is  found  not  to  be  of  great  value,  it  is  intended  to  omit  it  from  most  of 
the  computations  for  economy,  and  perhaps  later  to  restore  it  for  studying  its 
effects  on  the  best  filters  found.  If  doppler  is  omitted,  so  are  equations  (77}/ 
(78),  (79)#  (80)',  (87),  (88),  (89),  end  (90),  and  i  -  k-,5 ,6  in  (2^). 


3*5*4  Summary  of  the  computations.  The  sequence  of  computations  is  as  follows. 
First,  preliminary  investigations'  are  needed  to  establish  the  values  of  ov^, 

j  =  1,...,8,  in  (48)  and  of  t^,...,tw  for  each  trajectory. 


n 


I  1 

The  values  of  u^,  (44) ,  are  computed  for  the  reference  trajectory  and  the 
matrices  B^,  defined  in  (46)  and  computed  by  (48),  are  determined. and  stored. 


The  matrix  M  and  the  vector  d  are  computed,  by  (51)  through  (71),  for  the 
gun-radar  geometry.  M  and  d  are  used  in  (72)  through  (80) ,  together  with  the 
stored  values  of  u, ,  so  determine  the  projectile  state-vector  in  a  rectilinear, 1 
radar-centered  coordinate  system. 

t 

Then  (8l)  through  (100),  with  the  a-version  of.  equations  (8l),  (82),  (9l)j 
(92),  (93)#  (98),  (97)#  and  (98),  give  the  matrices  for  the  range -a- £  radar 

coordinate  system.  The  same  equations,  using  the  b-version,  gives  the  matrices 
for  the  range-azimuth-elevation  coordinate  system,  in  both  cases,  A^=  .• 

gives  the  matrics  needed  in  (26).  .  ■  .  1  : 

Then  for  the  spacing  tq  =  t^+1  -  t  of  radar  observations,  (26)  gives  the 

matrix  S;  (104)  is  used  for  any  other  evenly-spaced  radar  observations.  Finally 
(105)  through  (108)  give  the  covariance  matrix  P  of  the  least  squares  launch- 
point  estimator.  , 
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1 


l 


l 


I 
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b.d  "FILTER- SMOOTHER  DESIGN 


4.1  Introduction.  This  section  will  be  devoted  to  a  discussion  of  the  non- 
linear  filtering/smoothing  algorithm  for  estimating  instantaneous  state  vari¬ 
ables  that  will  be  implemented  by  the  University  in  the  next  quarter.  The 
first  subsection  'contains,  a  short  mathematical  description  of  the  algorithms 
to  be  employed.  The  second  wiil  desdribe  the  basic  problem  formulation.  This 
will  include  a  discussion  of  the  equations  of  motion,  state  variables,  and  co- 
variance!  matrices.  .Finally  the  third  subsection  will  describe _ implementation 

. concepts  and  future  plans.  ? 

'  '  .  '* 

•  •  *  .  i  .  . 

4.2  Estimation  Equations.  The  nonlinear  filter /smoother  specifies  estimates 

of  a  given  past  state  and  the  current  n-aimens^onal  state  of  a  dynamical 
system  observed  discretely  in  the  presence  of  additive  Gaussian  white  noise. 
The  system  is  assifmed  to  be  described  by  the  stochastic  vector  differential 
equation  ’  ,  , 


'  !  dxt  =  f(xt,t)dt  +  g(x.h,t)dp' 


V  '“Kt 


(109) 


where;  p,  represents  ai  Brownian  motion  process  with  covariance  matrix 
0  -  , 

■ 

I  .  *  , 

E(dpt,dp^}  =  ^  dt 


(no) 


The  corrupted  m-dimensional  (  ra  ^  n)  measurement  vectqr  z ^  is  related  to  the 
state  by 'the  expression  ■  : 


-  h<x^ + 


(UJ.) 


where 


Xx =  V  ,tmh- 


(112) 


V* 


(n3) 
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! 


I 


An  a  prion  estimate  of  the  state  x  with  covariance  cov(x  .x  )  is  also  assumed 
-  o  o  o 

given* 

The  algorithm  is  composed  of  a  set  of  discretized  extrapolation  equations 
and  a  set  of  estimate  update  equations.  The  extrapolation  equations  describe 
the  behavior  of  the  current  state  estimate  x3  and  the  associated  covariance 

Aj 

matrices  between  observation  samples.  The  update  equations  determine  improved 
estimates  of  -Jie  current  state  and  the  prescribed  previous  state  x^.(k  <  A). 

These  estimates  and  modifications  to  the  covariance  matrices  reflect  the  new 
information  available  in  z^. 

We  shall  denote  the  best  estimate  of  state  x^  based  on  a  realization  of 
measurements  Zj j  sfz^,.. .,z^}  as  x^^.  The  conditional  covariance  will  be-  de¬ 
noted  by  covCx^xJz^).  ■  '  '  -  ‘ 


The  extrapolation  equations  for 


*  ts-l  +  4t 

are  given  as  follows 


A*U  A 

00v(‘  ,XaIZA-1^  =  *£U-1  cov^Xj&-l,xA-llZje-l^Iu-l+  g^J6-m-l,tj8-l^A-l 

-  ‘  '  "  -  ^5) 

C°v(xk,xJZA_i)  .  cov(xk,x^1|ZjJ_1)5^_1  016) 

PA-lU-l  is  "the  nXm  Jac0^ian  matrix  of  the  f  (*  ,  t)  vector  function  with 
elements 

fo-Hje-J1'3  " 

3  *5 
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x 


^|4-1  =  I  +  At  F;_1|jUl  (XlB) 

Q^l  =  At  x)  019) 

where  I  is  the  nXn  identity  matrix. 

Defining  the  gain  term 

m  m  ”1 

=  c0V^x2,xjU (320 ) 

S  h. (x) 

3  =SX—  .  t321')' 

3  x  =  xm-i  ' 

We  have  for  the  update  equations 

4iW  ■  *»IM.  + 

\\t  *  \U-1  + 

°;>v<xW2e)  "  C0V<xJe>xilzi-i)(I  -  si  Hi  00',l-y-i’xi)h-i)) 

o°v(xk,xi|zJ1)  =  ewtv^lViX1  -  ^  h  cov<xi'xJ!lzi.l)) 

Air  expression  for  updating  the  approximate  covariance 
for  the  smoothed  estimate  also  exists.  The  expression  is  however,  not  required 
for  implementing  equation  123.  Its  utility  is  restricted  to  that  of  evaluating 
the  quality  of  the  estimate.  Since  algorithm  evaluation  will  be  an  important 
part  of  the  project,  we  include  the  equation 


(122) 

(323) 

(324) 
(125) 
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cov(xk,xklz£)  «=  c°v(xk,xk|ZjJ_1)  -  cov(xk,x£|Zji_1)Sjj  c^vV^-l^  (Og6) 


The  algorithm  is  initiated  by  the  given  a  priori  information 


A 


(227) 


and 


cot(x0,x0|Zo)  -  <*»(V*o> 


(128) 


4.3  Equations  of  Motion.  We  assume  ;he  radar  to  be  located  at  an  altitude  hQ  and 
at  laditude  n  on  a  rotating  spherical  earth  as  illustrated  in  Figure  2.  A  moving 


x  -  North 
g 

y  -  local  vertical 
g 


z  -  East 
g 


Fig  2  Radar  Location 


Cartesian  coordinate  system  x^y  ,z  will  be  defined  centered 

that  the  x  axis  points  north  end  the  y  axis  points  up  along 
S  g 


at  this  site  such 
the  local  verticle. 


Figure  3  illustrates  a  second  coordinate  system  xr>yr>zr 


Radar  faces  $  radius 
from  north 

Y  =  Tilt  angle 


Fig  3  Radar  Centered  Coordinate  Systems 


obtained  from  the  first  by  a  rotation  in  azimuth  and  y  in  elevation.  In 

this  system  the  x^  axis  is  normal  to  the  radar  face  pointing  outwards  and  the 

y  and  z  axes  are  in  the  plane  of  the  face, 
r  r 

The  simplified  dynamic  equations  of  motion  in  the  xr>yr>zr  coordinate 
system  are 


X  , 

r 

i 

>* 

1 _ 

vyr  gz  "  vzr  sy 

.  .  - 

X 

r 

V 

•  • 

yr 

=  -p  C  KB(Jf)v 

Vyr 

-S  K(*5) 

vzr  ^x  ”  vxr  sz 

+  205rU 

yr 

+ 

sy 

2 

•  • 

z 

V 

V 

v  g  -  v  „  g 

z 

g ,, 

r 

zr 

xr  &y  yr  ex 

- 

r 

drag  drift  coriolis  gravity 


(129) 
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and 


S  =  -.004  picA  (M)VS.  (130) 

where 

2  2  2  2  2 
h  =  h i  +  x .  siny  +  y  cosy  +(x  cos  y  +  y  sin  y  +  z  ) 

OX  X  X  X  X 

2(Ke+V 

H  =  Radius  of  the  earth 
e 

go  =  -9.80665  [1-0026  COS  2|i]  (Re/(Re  +  hQ))2 

=  g0[siny  +  {xr(l  -  3  Si  A )  f  yr(r'3‘  sistfcosy)}/^/*  hQ)]  \ 

gy  =  eo[cosy  +  {xr(-  3  sinVcosy  )  yr(l  -  3  cos2y)}/(Re  -<■  hQ)] 

gZ  =  go  2r/(He  +  V- 

The  skew  symmetric  corioulis  matrix  is  defined  by  the  upper  off-diagonal 

terms 


UU(J  =  0  sin§  cosji 

0 


cos$sinycosn  -  cosysinji 

cosSPcosycosfi  +  sicysiny. 
0 


(131) 


and  the  projetile  velocity  relative  to  the  air  is  given  by  the  vector  (of 
magnitude  V): 
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V  =  x  -  W  (cos  AZ  cosy  cos§  +  sin  AZ  cosy  sin§) 
xx*  r  s 


V  =  y  -  W  (-cos  AZ  sincy  cos§  -  sin  AZ  siny  sin?1) 
yr  r  s  ' 


Vzr  =  \  ”  Wg(cos  AZ  sin?  -  sin  AZ  cos§)  (132) 


where  Wg  is  the  magnitude  of  wind  and  AZ  its  azimuth  direction  (measured 
clockwise  from  north). 

Atmospheric  data  are  assumed  available  in  the  look-up  format  (as  a  function 
of  h)  prescribed  for  the  BRL  pointraass  program.  The  universal  aerodynamic 
curves  are  described  by  piecewise  fourth  order  polynomials  of  mach  number* 


^(M)  =2^ 

K(M)  =  2  b.  M1  (133) 

ka(m)  =  s  Si  m3- 


4.4  State  Variables.  As  indicated  in  the  previous  quarterly  report,  the  un¬ 
known  aerodynamic  parameters  will  be  treated  as  additional  state  variables  to 
be  identified  by  the  estimation  equations  122  and  123.  Included  in  the  simpli¬ 
fied  dynamics  (eq,  129)  are  two  such  parameters  C  and  S.  The  variable  S  is  the 
product  of  two  unknowns,  an  assumed  locally  constant  drift  parameter  r  similar 
to  C  and  the  projectile  spin  rate. 

The  state  vector  x^  is  defined  by  the  expressions 


*4 


*7 

x8 


(134) 


And  from  the  equations  of  motion  we  obtain 


f(x)  = 


DlVx 

D1Ty 

D1V 

z 


*  K<Vz 

+  ffi(vA 


VZSy>  +  2a>K(U12:S+U13X6)+^ 

\%>+2“r'U21X4+U23  V+gy 

Vy^  +  2f“R(U«X4+U32X5)+^ 

0 

D3  x8 


(135) 


where 

D1  =-p  ic(M)V 

L2  =  -  X8  K^ltf/V2  C^6) 

D3  =  -  -004p  Ka(M)V 
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The  components  of  the  J'acobian  matrix  F  are  given  in  Appendix  E. 


The  x^,  z^  coordinate  system  was  orginally  selected  because  the  phased 

array  radar  measurements  are  defined  directly  in  this  system  as  is  illustrated 
in  Fig.  4. 


Fig  4  Observation  Vector 


In  state  variable  notation  these  measurements  (including  optional  doppler  rate 
information)  are 


/R  ‘ 

r  2  2  2  1/2 

(xx  +  x2d  +  x3") 

h(x)  = 

a 

= 

cos"1(X1/R) 

(3 

cos_1(X3/R) 

R 

(Xi  X4  +  X2  X5  +  X3  X8)/Ra 

u 

Again  the  Jacobian  matrix  components  are  given  in  Appendix  E. 
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4.5  Covariance  Terns.  In  radar  tracking  applications  the  Brownian  motion 
process  in  equation  1  is  used  to  represent  random  atmospheric  effects  and 

errors  due  to  modelling  approximations.  Of  the  latter,  the  most  significant 
is  the  assumption  that  C  and  r  are  constant  functions  of  mach  number.  This  is 
particularly  true  for  the  drift  constant,  however  as  illustrated  in  Pig  5  a 
constant  C  is  also  not  an  especially  good  approximation  for  K^_.  (M)/k^(m)  in 

the  region  about  mach  1. 

Although  it  has  been  shown  that  a  best  Q  matrix. for  filter  implementation 
must  be  determined  experimentally,  order  of  magnitude  estimates  for  and 

may  be  determined  as  follows.  We  assume  that  at  any  time  t  for  particular 

shell  j  C  has  been  estimated  such  that 


C  =  Kp  (mJ/k^M). 


(138) 


Then  at  time  t+At  due  to  the  change  in  mach  number  (and  thus  the  true  value 
of  the  ratio  of  equation  (138)  there  will  be  an  error  introduced  in  C 

given  by  the  expression 


‘JD  “  It  »D 

■MskwW 

Vs  3 

“At  TTQm)  [  3m  *b  cm)  -  ca  k;<m)] 

Vs  D  d  “35J — 

Recalling  the  assumption 

K  M  _  s  a3 

D1  0,4  1 


(139) 

(140) 

(141) 

(142) 


we  obtain 
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(143) 


v 


L. 


Da.).  M1"1 

12. 


where  we  approximate 


V  =  jp  Kjj  (M)  V  V  +  gj 


(144) 


Rather  obviously  the  development  for  the  drift  coefficient  will  be  very  similar 

to  the  above.  By  computing  (for  nominal  p,  V  ,  and  QJ3)  error  terms  for  drag 

s 

and  drift  as  a  function  of  mach  number  for  the  given  standard  shells  and  then 
averaging  the  results;  the  following  tentative  g(x)  8x5  and  Q‘5X5  matrices  were 
determined 


r 


.0016 


.0016 


.0016 

%k 

%5 

_ 


1 


D3X8 


(145) 


(146) 


where 


044 


j- (13.634  -  12.34m) 
.-(29MU  +-28.58) 


1.14  <  M  ;  M  £  .7 


.7  <  M  ^  1.01 
1.01  <  M  £  1.14 
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000196 


1.10  <  H  }  M  ^  .83 


%5=  1 o-03-382  -  n.666M)  .83  <  M  i  .96 


(147) 


.  -  (lOM  -  7.2 


.96  <  M  ^  1.10 


4.6  Implementation.  The  filter/smoother  will  be  simulated  using  data  gener¬ 
ated  by  the  programs  described  above  in  Section  2  of  this  report  for  a  number 
of  nominal  projectile  trajectories.  The  past  state  to  be  estimated  by  fixed 
point  smoothing  will  be  that  which  would  correspond  to  Xq,  the  first  point  of 

the  tracked  portion  of  the  trajectory.  Experiments  will  be  conducted  as  to 
the  viability  of  using  algorithm  generated  covariance  matrices  for  error  evalua¬ 
tion  studies.  Convergence  properties  of  the  algorithms  relative  to  given  Q 
matrices  and  higher  order  expansions  (particularly  h(x))  will  be  considered. 
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5.0  CONCLUSIONS  AND  PLANS  FOR  FUTURE  WORK.  As  pointed  out  in  Section  2  above, 
a  set  of  simplified-dynamics  equations  has  been  obtained  and  programmed  with' 
corrections  for  both  drag  and  drift.  Tentative  numerical  results  obtained  so 
far  indicate  that  these  equations  approximate  rather  well  the  effects  of  drift 
keeping  the  resultant  lateral  errors  quite  small.  The  drag  term,  however  leads 
to  longitudinal  errors  that  are  disappointingly  large.  This  effect  is  most  pro¬ 
nounced  for  cases  where  the  point  at  which  backtracking  begins  occurs,  after  the 
projectile's  velocity  has  crossed  the  mach  1  region.  .This  appears  to  indicate 
that  satisfactory  results  can  only  be  obtained  in  such  cases  when  it  is  possible 
to  obtain  valid  tracks  of  the  target  while  it  is  still  going  at  speeds  greater 
than  mach  1.  Work  remaining  to  complete  this  effort  must  include  a  systematic 
investigation  of  the  lateral  and  longitudinal  errors  to  be  .expected  in  backtrack* 
ing  for  a  comprehensive  set  of  trajectories.  Also,  additional  thought  will  be 
applied  to  possible  refinement  of  the  drag  correction  so  as  to  decrease  the  long-  , 
itudinal  errors  to  an  acceptable  value.  '  .  :  ! 

1 

In  Section  3  the  formulation  necessary  to  determine  the  magnitudes  of  irre- ' 
ducible  errors  due  to  radar  observation  noise  has  been  developed  and  the  necessary 
programming  is  underway.  The  purpose  of  this  progrm  is  to  determine  the  theoreti¬ 
cal  accuracies  attainable  in  launch-point  estimation,  on  the  assumptions  that  the 
simplified-dynamics  equations  are  correct.  That  is, -it  finds  the  minimal,  errors, 
using  best  possible  filtering  and  smoothing  without  limits  on  the  amount  of  qompu- . 
tution  required,  that  result  from  the  existence  of  the  radar  noise.  These  minimal 
errors  resulting  from  radar  noise  of  course  add  to  the  errors  introduced  by  use  of 
the  simplified-dynamics  equations.  This  will  include  a  determination  of  whether 
acceptable  estimates  can  be  made  for  the  drift  and  drag  parameters  so  as  to  permit 
effective  use  of  the  simplified-dynamics  equations  for  backtracking.  .  1 

i  * 

The  work  in  Section  4  constitutes  a  detailed  formal  procedure  for  performing 
the  filtering  and  smoothing  required  to  estimate  necessary  inputs  for  backtracking 
from  a  set  of  radar  observations.  Work  has  begun  to  create  a  computer  program  for 
the  realization  and  evaluation  of  this  algorithm.  ’  • 
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Appendix  A 


Modified  Point-Mass  Equations 


The  equations  here  described  are  taken  from  Reference  [l] .  The  total 
instaneous  vector  acceleration  of  the  projectile’s  center  of  mass  with  re¬ 
spect  to  the  ground  is 


u  =  (drag  term)  +  (drift  term)  +  (Magnus  term) 

+  (gravity  term)  +  (Coriolis  term). 


The  drag  uerm  is 


■(pd^/m)^  +  Kp  (Q0e)2}  v  v  , 


the  drift  term  is 


(pd^ /m)KjV\£  a£ 


and  the  Magnus-effect  term  is 
(pd3/m)KpNQ(de  x  v) 


where  p  is  the-  air  density, 

d  is  the  projectile  diameter, 
m  is  the  projectile  mass, 

Kjj  ,  Kp  ,  K^,  and  K^,  are  functions  of  Mach  number  characteristic 
of  the  projectile, 


Q  is  the  projectile  yaw  drag  factor, 

*  l  is  the  projectile  lift  factor, 

N  is  the  spin  angular  velocity, 

v  is  the  vector  of  projectile  velocity  with  respect  to  air, 
and  a  is  an  approximation  for  the  angle  of  repose. 

— ♦  c 

Auxiliary  equations  are  used  for  computing  N  and  a  : 
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N  =  (pd  /A)K.Nv 


where  A  is  the  projectile's  axial  moment  of  inertia  and  K.  is  a  function  of 
Mach  number.  A 

Se*  (ab  -aa)(^xH)  ~%(ZXS) 

where  g  is  the  gravity  vector, 

aa  =  Kjj  v4  +  pd5  Kp  Kt  K2/) 

0<b  =  mi^N/CpdK^  Kj^  v4  +  pd3  Kp  K^J,  mV  ) 

and  Ky,  are  functions  of  Mach  number. 

From  data  not  given  in  Reference  [l]  but  used  in  the  programs  that  com¬ 
pute  the  modified  point-mass  equations  for  the  105mm,  155mm,  175mm  and  8 
inch  projectiles,  =  0,  giving  at  once  the  simplification 

«e  =  "(AN/pd^Kv  x  u) 


and  the  drift  acceleration  term  becomes 


- (An/mdK^v2 )KjA (v  x  u). 

In  2.2,  is  called  K^(m)  and  AA/md  is  called  a.^. 

Computation  of  the  urag  coefficients  and  K^Qal)  indicate~that  for 
ae  £  2°  the  sum  behaves  to  an  order  of  magnitude  as  *  alone.  Since  Simula- 

DO 

tions  with  the  modified  point  mass  model  indicate  that  angles  of >2  occur  only 
in  the  vicinity  of  apogee  wad  then  only  for  QE's  £  800  mils,  it  is  assumed  that 
the  approximation 

H  =  (-p/C)!^  v  % 
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is  always  valid. 


Appendix  B 


Universal  Aerodynamic  Functions 

This  appendix  includes  coefficients  for  and  plots  of  the  universal  drag  and 
drift  functions  employed  in  the  estimation  algorithms.  For  comparison;  graphs 
of  'the  nominal  drag;  and  drift  aJK^M)  functions  for  the  four  standard 

projectile  types  employed  in  our  analysis  are  also  included.  It-  can  easily  be 
seen  that  the  universal  d^ift  function  tends  to  be  a  poorer  approximation  than 
does  the  corresponding  drag  function. 

The  multiplying  constants  c..  and  determined  with  the  universal  curves  are 
given  along  with  their  respective  nominal  ballistic  coefficients  and  lift 
factors  &.  in  Table  3.1.  Computation  of  the  appropriate  and  r^  for  any  part¬ 
icular  trajectory  (charge  and  Q.E. )  in  the  backtracking  evaluation  required  the 
x'ollowing  scaling  .operations 

ci  -  °i  (ci  /°i  > 

true  nom  nom  true 

and 

r.  =  r.  (i.  /&.  ) 

1  i  '  i  ’  x 

‘true  nom  true  nom 


where  the  true  and  JL  true  are  available  from  BRL  data. 

Table  B.l 


Shell 

Diameter 

Drag 

1 

Nominal 

Ballistic  Coef. 

Drift 

r. 

1 

Nominal 
Lift  Factor 

'  105  MM 

1.38850 

1.919 

.661707 

.863 

155  MM 

1.07657 

2.331 

1, 12855 

.963 

175  MM 

.69031 

3.101 

.7411013 

1.009 

8  IN 

.84459 

3.16 

1.468419 

.880 

The  coefficients  for  all  three  aerodynamic  functions  also  given  in  this  appen 
dix  as  Table  B.2  are  for  piecewise  fourth  order  polynomials  (of  Niach  number)  fits 
to  the  respective  curves.  The  particular  coefficients  are  assumed  valid  between 
the  previous  break  point  mach  number  up  to  ar.d  including  the  value  given  in  column 
six  of  the  particular  row  of  coefficient  values.  All  values  are  given  in  the  same 
format  as  is  employed  in  the  "aeropacks 11  of  the  BKL  point  mass  program. 

last  of  Figures 

Fig  B.l  Universal  Brag  Curves 
Fig  B.2  Zero-Yaw  Drag  Curves 
Fig  B.3  Universal  Drift  Curve 
Fig  B.4  Zero- Yaw  Drift  Curves 
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Apnendix  C 


Computer  Program  Development 

The  starting  point  for  this  work  was  a  FORTRAN  program  oo  compute  modified- 
point-mass  trajectories,  produced  at  Ballistics  Research  Laboratory  and  furn¬ 
ished  to  Project  RATRAN  by  USAECOM.  In  adapting  this  program  to  the  Spectra-70 
Computer  at  the  Moore  School  several  minor  bugs  were  eliminated  and  several 
features  were  added  to  obtain  desired  results.  A  description  of  the  main  pro¬ 
grams  in  use  is  as  follows. 


F-l.  PTMASS  -  a  module  composed  of  9  subroutines  which  computes  the  trajectory 
of  a  shell  and  the  corresponding  radar  coordinates  in  either  phased-array 
or  range-azimuth-elevation  form.  This  module  includes  the  following  sub¬ 
routines  : 

A.  PTMASS  -  the  main  calling  program  which  calls  all  of  the  others  and  which 
reads  in  and  stores  the  variable  data  -  corresponding  to  a  particular  tra¬ 
jectory.  This  input  data  includes  the  necessary  'variables  for  computation 
of  the  radar  coordinates. 

B.  AEP.OBL  -  a  subroutine  called  by  PTMASS  which  reads  the  aeroballistic  data 
for  a  particular  shell  from  a  disk  file,  checks  the  data  for  errors,  and 
stores  it  for  later  use  in  computation  of  the  aerodynamic  coefficients  KDA, 
KDO,  KLO,  KM,  KF,  KA,  and  KLA. 

C.  TRAJ  -  a  subroutine  called  by  PTMASS  which  computes  the  velocity  and  posi¬ 
tion  vectors  of  the  trajectory  at  specific  points  in  time,  using  a 'single- 
step  predictor-corrector  type  of  numerical  integration,  and  which  also 
computes  the  corresponding  radar  coordinates.  It  is  possible  to  obtain 
either  phased-array  or  range-azirauth-elevation  radar  coordinates  by  vary¬ 
ing  an  input  variable.  This  subroutine  will  compute  the  trajectory  until 
one  of  the  various  "stop  conditions"  is  satisfied.  These  stop  conditions 
include  termination  of  computation  when  time,  the  x-coordinate  of  the 
position  vector,  the  range,  the  height  of  the  up- leg,  or  the  height  of 
the  down-leg  have  reached  a  given  value,  or  at  the  summit  of  the  trajec¬ 
tory.  Other  values  printed  out  at  each  time  interval  include  mach  number, 
drag,  lift,  and  magnus  force. 

D.  COMPC  -  a  subroutine  called  by  PTMASS  before  TRAJ  is  called  to  compute  the 
trajectory.  COMPC  computes  the  ballistic  coefficient  c  as  a  function  of 
the  quadrant  elevation. 

E.  COMPT  -  a  subroutine  called  by  PTMASS  and  TRAJ  which  computes  the  time  of 
flight  correction  as  a  function  of  machine  time. 

F.  COMPFS  -  a  subroutine  called  by  PTMASS  which  computes  the  fuze  setting. 


G.  COMPL  -  a  subroutine  called  by  PTMASS  before  the  trajectory  is  computed 
which  computes  the  lift  factor  as  a  function  of  the  quadrant  elevation. 
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H.  COMPNF  -  a  subroutine  called  by  PTMASS  which  computes  the  muzzle  velocity 
correction  factor  N  if  the  standard  weight  of  the  projectile  is  different 
from  the  actual  weight. 

I.  COMPPT  -  a  subroutine  called  by  PTMASS  which  computes  an  adjustment  to 
the  initial  velocity  if  propellant  temperature  is  not  70°. 

F-2.  BACK. PTMASS  -  a  module  of  10  subroutines  which  computes  the  trajectory  of 
a  shell  and  has  the  option  of  stopping  at  any  point  in  time  and  integrat¬ 
ing  backwards  until  height  is  zero  using  the  simplified  dynamics  equations 
for  drag  and  drift.  This  module  includes  many  of  the  same  routines  as 
PTMASS  -  those  that  are  different  are  explained  in  further  detail  below. 

A.  AGAIN  -  the  main  executive  which  calls  all  the  others.  AGAIN  is  Similar 
to  the  routine  PTMASS  except  that  it  does  not  have  an  insert  allowing  it 
to  read  radar  data.  It  also  differs  from  PTMASS  in  that  it  has  a  patch 
which  calls  AEROBL  to  read  in  the  universal  aeroballistic  pack  for  use 
in  backward  integration. 

B.  AEROBL  -  same  as  before 

C.  AGAIN. TRAJ  -  a  subroutine  called  by  PTMASS  to  compute  the  trajectory.  This 
routine  is  similar  to  TRAJ  except  that  it  does  not  compute  corresponding 
radar  coordinates  at  each  point  of  the  trajectory. 

D.  BTRAJ  -  a  subroutine  called  by  PTMASS  which  computes  the  trajectory  back¬ 
wards  until  height  is  zero  using  the  simplified  dynamics  equations  for 
drag  and  drift. 

E.  COMPC  -  same  as  before 

F.  COMPT  -  same  as  before 

G.  COMPFS  -  same  as  before 

II.  COMPL  -  same  as  before 

I.  COMPNF  -  same  as  before 

J.  COMPPT  -  same  as  before 

F-3*  AGAIN. PTMASS  -  a  module  of  10  subroutines  which  computes  the  trajectory 
of  a  shell  until  one  of  the  various  stop  conditions  is  satisfied.  This 
program  has  the  option  of  stopping  at  any  point  in  time  and  integrating 
backwards  until  height  is  zero  using  the  standard  BRL  equations  for  drag, 
drift,  and  magnus  force.  This  capability  was  designed  into  the  module  as 
a  check  on  the  backward  computations  of  BACK. PTMASS  with  the  simplified 
dynamics  equations. 
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A.  AGAIN  -  same  as  before 


3.  AEROBL  -  same  as  before 

C.  AGAIN .TRA  J  -  same  as  before 

D.  AGAIN. BTRAJ  -  a  subroutine  called  by  PTMASS  which  computes  the  trajectory 
backwards  until  height  is  zero  using  the  standard  equations  for  drag,  drift, 
and  magnus  force. 

S.  COMPC  -  same  as  before 

F.  COMET  -  same  as  before 

G.  COMPFS  -  same  as  before 

H.  COMPL  -  same  as  before 

I.  COMPN?  -  same  as  before 

J.  COMPPT  -  same  as  before 


F-U.  Considerations  in  implementing  simplified  drag  and  drift: 

As  originally  designed,  the  program  module  AGAIN .PTMASS  had  the  capability 
of  stopping  computation  of  a  trajectory  at  any  point  in  time  [using  a  stop  code 
of  U]  and  of  integrating  backwards  until  height  was  zero  [using  a  stop  code  of 
9] •  It  wo> ild  have  been  easy  to  run  the  program  with  a  stop  code  of  4  to  some 
point  in  time  and  then  re-run  the  program  using  the  final  values  of  the  velocity 
and  position  vectors  from  the  previous  run  as  input  to  compute  backwards  with  a 
stop  code  of  9.  This  approach,  however,  would  have  been  wasteful  of  computer 
time  and  slower,  as  the  entire  program  would  have  to  be  reloaded  to  compute  back¬ 
wards,  and  the  input  values  for  the  second  half  of  the  run  would  have  had  to  be 
written  to  a  disk  file,  or  punched  into  cards.  Instead  the  module  was  modified 
so  that  upon  raading  a  certain  input  card,  the  program  automatically  switches  to 
a  subroutine  which  integrates  backwards  using  the  simplified  dynamics  equations, 
taking  the  last  values  of  the  trajectory  as  a  starting  point. 

The  procedure  is  initiated  by  reading  an  A  card  in  the  input  data.  This 
signals  the  program  to  read  in  the  universal  aeroballistic  pack  for  that 
particular  shell  type.  Then  the  subroutine  BTRAJ  is  called  to  perform  the 
backward  integration.  The  patch  made  to  the  program  is  shown  below  in  lines 

620-633. 
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oo  o  u  o  uooo 


GOTO  110 

PROCEDURE  FOR  AN  A  CARD 
STORE  NEW  AEROSALLI STIC  OATA 

112  IDSET=B( 1 ) 

TEST  READ  NEW  AEROPACK  OR  NOT 
IFtIDSET.EQ.32)  GO  TO  700 
CALL  AEROBL(IDSET) 

PRINT  HEADINGS 

700  IF(TEST.EQ.O)  GO  TO  701 
WRITE(6»406) 

GO  TO  702 

701  WRITEt  6*405 ) 

702  HRITE(6*450) 

COMPUTE  THE  BACKWARD  TRAJECTORY  WITH  SIMPLIFIED  EQUAT 
CALL  BTRAjtl, METRO) 

GO  TO  (43*113), I ER 

113  GOTO (95* 109, 108 ) *  I ERR 


The  subroutine  BTRAJ  was  created  from  the  rou+  TRAJ  which  computes 
trajectories  under  a  variety  of  stop  conditions  A  patch  was  made  just  after 
the  entry  point  to  the  routine  which  sets  all  the  necessary  variables  to  en¬ 
able  backward  integration  to  take  place  (stop  code  =  9)*  It  was  necessary  to 
remove  all  statements  which  initialize  the  values  of  time,  position,  and 
acceleration  for  the  final  values  from  the  forward  run  were  to  be  used.  Since 
these  variables  were  all  placed  in  the  COMMON  area  of  storage,  the  final  values 
were  automatically  passed  in  when  the  subroutine  was  entered.  The  patch  that 
was  added  to  set  the  variables  is  shown  below  in  lines  44-52  of  the  listing* 
DATA  TOL/  1.01/  ~  - 

METR®METRO 
GOTO ( 1000*2 ) , JTR 

C  FIRST  CALL  PROCEDURE 

C  SET  TIME  TO  FINAL  VALUE  FROM  FORWARD  RUN 

1000  Yl=T 

C  SET  STOP  CODE  TO  9 

IST0P=9 


C 


SET  PRINT  INTERVAL  TO  -1.0 


PINT®- I.  -  . 

C  SET  FINAL  VALUE  TO  ZERO 

FV®0. 

RECW=2 .20462 2622 /WT 

Other  changes  made  to  the  program  to  implement  simplified  drag  and  drift 
include  a  section  to  set  up  the  constants  used  in  the  actual  computation  of 
drag  and  drift,  as  explained  in  another  section  of  this  report,  and  replace¬ 
ment  of  the  equations  to  compute  the  forces  of  drag,  drift,  and  raagnus  force 
with  the  new  equations  as  shown  below: 


ALPHA® SORT ( ALPHS ) 

FIND  COMPONENTS  OF  UDOT— THE  ACCELERATION 
D1  IS  DRAG 
02  IS  ORIFT 
D3‘  IS  MAGNUS' 'FORCE 
710  D1=-RH0*V#KD0*DDD 
D2  ®  NR*KL0/VSQ 
03=0* 
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Appendix  p 


Coordinate  Conversions  '  '• 

Project  RATRAli  has  been  using  a  computer  program  written  by  ERL,  for  generat¬ 
ing  projectile  trajectories  by  means  of  the  modified  point-mass  equations.  In 
order  to  prepare  for  simulating  the  signal-processing  filters  to  be  designed, 
there  have  been  added  instructions  to  the  BRL  program  that  compute  the  eenter-of- 
mass  projectile  coordinates  in  radar  axes,  concurrently  with  the  computations  by 
the  original  program  of  the  position  in  gun-axes.  This  appendix  describes  the 
radar-axas  computations.  ’ 

The  gun-axes  in  the  original  program  are  a  right-handed  rectilinear  system 
with  origin  at  the  gun,  the  y-axis.  positive  upward  along  the  local  vertical,  aud 
the  initial  velocity  vector  in  the  x-y  plane.  The  radar  coordinates  ol  the  center 
of  mass  of  the  projectile  are  computed  by  rigid  translations,  and  rotations  of.  cor 
ordinate  systems  until  the  projectile  position  is  expressed  in  a  right-handed  ; 
rectilinear  system  with  the  origin  at  the  radar,  the  y-axis  positive  upward  along 
the  local  vertical,  and  the  x-axis  horizontal  along  the  nominal  zero-azimuth  line 
of  the  radar.  Then  the  rectilinear  coordinates  are  transformed  to  range -a- $  radar 
coordinates  and  also  to  'range-azimuth-elevation.  The  earth  is  assumed  to  be  a 
sphere  of  radius  R.  The  radar  elevation  is  not  assumed  to  'be  the  same  as  the  gun 
elevation.  ■  1 

Three  input  constants,  xg,  yg,  and  A2,  define  the  gun-axis  system  with1 respect 

to  the  radar  coordinates.  These  three  constants  are  most  easily  explained  on  the 
basis  of  a  flat  earth.  A2  is  measured  clockwise  looking  down,  from  a  line  through1 
the  gun  parallel  to  the  radar  x-axis,  as  shown  in  figure  D-l.  The  position  of  the 
gun  in  the  radar  rectilinear  coordinate  system  is  at  z  =  x  ,  x  =  y  .  {The  choice 

of  the  symbols  xg,  yg  here  and  in  the  program  annotations  fs  inconsistent  with  the 

choice  of  axes,  for  historical  reasons  arising  in  the  fact  that  information  from 
BRL,  including  computer  programs,  uses  a  vertical  y-axis  and  information  from  ECOM, 
in  particular  the  gun  location,  uses  a  horizontal,  y-axis. )  ; 

The  first  coordinate  transformation  is  a  rigid  rotation  about  the  vertical,  to 
make  the  gun  x-y  plane  pass  through  the  radar.  The  equations  are 

i  ' 


Cl  =  (-x  sin  A2  -  y  cos  A2)/d  ’ 

®  s 

Si  =  (  x  cos  A2  -  y  sin  A2)/d  , 

8  8  > 

i 

i 


,  2  2,1/2  ’ 
(xg  +  ye 


I 


I 


t 


I 


x..  =  x  Cl  -  z  SI 

1  O  ■  O 


=  x  SI  +  z  Cl 
1  .o  o  - 


where  y.  ,  yQ,,  z  are  the  coordinates  of  the  projectile  center  of  mass  as  produced 
by  the  original  BRL  computer  program. 

!  I  i 

Next  there  is  a  rigid  rotation  about  the  center  of  the  earth,  to  make  the  y- 
axis  pass  through  the  radar;  ^he^origin  is  also  moved  upwferd  a  distance  h,  the 
height  of  the  radar  above  the  elevation  of  the  gun.  The  equations  are 


■'  e  =  q/r, 


!x2  =  Xx(i  -  e2/2)  -  yoe  -  d 


C3  =  -y 


=  x2C3  -  2^3  ! 

.  .  z3  =  x2S3  +  z1Cs  ( 

■!  » 

1 

^3  =  xl0  +  yo(1  "  02/2)  "  ,6d/2  "  h> 


where  h  is  the  height  of  the  radar  above  the  altitude  of  the  gun.  The  coordinates 
of  the'  ceriter  Qf  mass  of  the  projectile  ip  the  rectilinear  radar-axe’s  system  are 
i  X3,  y^,  *  .  The  reason  for  hot  explicitly  performing  the  rotation  through  0  in  the 

usual  manner  is  to  avoid  the  round-off  error  of  subtracting  nearly-equal  quantities 

of- the  ordhr  of  R,  the  radius  of  the  earth. 

■  ,  ■  :  1  .  '  1 

Concurrently  with  the  computation  of  x_,  y0>  z  there  are  computed  their  time 

'  l  J  #  ^ 

derivatives,  using  the  time  derivativex  x  ,  y  ,  and  z  produced  by  the  original 

I  *  OOO 


BRL  program: 


x,  =  x  Cl  -  z  SI 
1  o  o 

z,  =  x  SI  +  z  Cl 
1  o  o 

x2  =  ixa  -  e2/2)  -  yoe 

=  XgC3  -  3^3 
z3  =  x2S3  +  zxC3 

The  radar  range 

/  2  2  2,1/2 
r  =  (Xj  +  yj  +  ) 

and  the  range  rate 

f  =  (x3  x3  +  y3  y3  +  z3  z3)/r 

are  computed. 

The  phases-array  radar  angles  are  computed  by 

a  =  arc  cos(x3/r) 

0  =  arc  cos(z3/r) 

The  azimuth  and  elevation  angles  are 

AZ  =  arc  tan(z3/x3) 

EL  =  arc  sin(y3/r) 
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Radar 


'3 


Position  of  Gun  in  Rectilinear 
Radar  Axes 


(Note:  the  gun  coordinates  are  called  (x  ,y  )  as  shown  above,  in  the  program 

8  £ 

described  in  Appendix  D.  They  are  called  (z  ,x  )  in  section  2. ) 

S  8 


Figure  D-l 


Appendix  E 

Components  of  Jacobian  Matrices 
A.  F  Matrix 
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p45  =  A  ST +  sr  <Y* 
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f58  *  B2<Vx  -  WA) 
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-°(h)  x7  Lo>(1+1)  ai 
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Assume  U.S.  Standard  Atmosphere 


B.  H  Matrix 
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Hn  =  Xx/R 

H12  =  X2/R 
By  -  X3/r 


Hgi  =  - 


^2  = 


R2  Ag2  +  X  2 


X1  X3 

^3  =  „2  X,  2  .  „  2 


R“  -/Xg  +  X3£ 


X3  X1 


*»‘77g7^ 


X  X2 

**32  3  2  /  2  2 
J  R  /X^  +  Xg 


i 

I 

* 

i 
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Appendix  F 
Glossary  of  Symbols 


0 f,p 
AZ 
C 
d 

Dl 


Ifi 

D3 

f(xt,t) 

FA-1|A-1 

So 

*x>VSz 

g(x,t) 

Y 

h 

h 

o 

h(Xjj) 

Kj(M) 

\(M) 

k^m) 


Projectile  Yaw  angle  of  repose 
Radar  angle  measurements 
Wind  azimuth  angle 
Unknown  drag  parameter 
Projectile  diameter 


-  p  C  K(M)V 

-  S  K^Mj/V2 

-  .004  p  k*(m)v 

State  variable  equations  of  motion 

Jacobian  matrix  df(x)/dx 

Sea  level  acceleration  of  (gravity 

Gravity  vector  components 

Matrix  function 
Radar  tilt  angle 

Target  altitude 
Radar  site  altitude 
Radar  measurement  vector 
Jacobian  matrix  dh(x)/dx 
Drift  function" 

Drag' function  C  Projectile  i 


Spin  function 
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Universal  functions 


K(M),  KD(M),  Ka(M) 
1 

m 

M 

P 

N 

®R 

i 


% 


r 

P 

R 

t 

R 

ft 

s 

u 


Projectile  lift  factor 
Projectile  mass 
Mach  no.=  v/Vs 
Laditude  of  radar 
Projectile  spin 
Earth  rotational  rate 

Radar  pointing  angle  (clockwise  from  north) 
1  +  At  F£-l|£-i 

State  noise  covariance  matrix 

V* 

Unknown  drift  parameter 

Density  of  air 
Radius  of  earth 

Target  range 

Range  rate 

NXr 

Coriolis  matrix 
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Projectile  velocity  (W.R.T.  air) 
projectile  velocity  vector  components 
Speed  of  sound 

Ground  located  Cartesian  coordinate  system 
Radar  located  cartesian  coordinate  system 
State  vector 
Estimated  state  vector 


Radar  measurement  vector 


